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Methods of quantum nuclear wave-function dynamics have become very efficient in simulating large isolated 
systems using the time-dependent variational principle (TDVP). However, a straightforward extension of 
the TDVP to the density matrix framework gives rise to methods that do not conserve the energy in the 
isolated system limit and the total system population for open systems where only energy exchange with 
environment is allowed. These problems arise when the system density is in a mixed state and is simulated 
using an incomplete basis. Thus, the basis set incompleteness, which is inevitable in practical calculations, 
creates artificial channels for energy and population dissipation. To overcome this unphysical behavior, 
we have introduced a constrained Lagrangian formulation of TDVP applied to a non-stochastic open system 
Schrodinger equation (NOSSE) [L. Joubert-Doriol, I. G. Ryabinkin, and A. F. Izmaylov, J. Chem. Phys. 141 , 
234112 (2014)]. While our formulation can be applied to any variational ansatz for the system density matrix, 
derivation of working equations and numerical assessment are done within the variational multiconfiguration 
Gaussian approach for a two-dimensional linear vibronic coupling model system interacting with a harmonic 
bath. 


I. INTRODUCTION 

Modelling photo-induced quantum dynamics in large 
molecular systems requires consideration of relatively lo¬ 
calized chromophoric regions where this dynamics orig¬ 
inates as well as macroscopic environment embedding 
these regions and influencing their behavior. 1-3 This is a 
very challenging task due to exponential scaling of quan¬ 
tum mechanics with the number of degrees of freedom 
(DOF). 

One successful approach to reduce the prefactor of the 
exponential scaling for isolated systems in pure states 
is to project the time-dependent Schrodinger equation 
(TDSE) onto a time-dependent basis set. The basis time- 
dependence gives flexibility to minimize an error in an 
approximate solution of the TDSE with reduced compu¬ 
tational cost as compared to that for static bases. Equa¬ 
tions of motion (EOM) for the basis set parameters are 
obtained using the time-dependent variational principle 
(TDVP). 4-6 This approach led to methods that are able 
to describe up to hundreds of nuclear DOF such as multi¬ 
configuration time-dependent Hartree (MCTDH) single- 
7,8 and multi-layer 9 methods, as well as variational mul¬ 
ticonfiguration Gaussian (vMCG) method. 10,11 However, 
computational costs of these approaches still scale ex¬ 
ponentially with the number of DOF so that molecular 
systems such as proteins (several thousands of DOF) can¬ 
not be simulated yet. Moreover, wave-function methods 
cannot describe energy dissipation and quantum decoher¬ 
ence of a system interacting with environment in a mixed 
state (e.g., heat reservoir or incoherent sun light 12 ) with¬ 
out departing from the TDSE. 13 

The quantum master equation (QME) formalism in¬ 
volving the system density matrix (DM) can model in¬ 
teraction of the quantum system with large macroscopic 
environment and describe associated processes of energy 


exchange and quantum decoherence. 13,14 These processes 
give rise to mixed states of the quantum system, which 
cannot be described with a single wave-function. There¬ 
fore, the use of a DM, /5(f), becomes indispensable, and 
the TDSE is replaced by the QME 

p(t) = C[p(t)\, (1) 

where the dot stands for the partial time derivative, and 
£ is the Liouvillian super-operator. In the case of an 
isolated system with the Hamiltonian H, 

£[p(t)] = -i[H,P{t)], (2) 

and Eq. (1) becomes the Liouville-von Neumann equation 
that generates a unitary evolution of the DM in which 
the system energy is conserved. For an open system, £ 
also contains a non-unitary part that is responsible for 
energy dissipation and quantum decoherence. This non- 
unitary part is obtained either from purely mathematical 
reasoning to preserve semi-positivity of the system den¬ 
sity (e.g., Lindblad approach 15 ) or derived using pertur¬ 
bation theory with respect to the system-bath coupling, 
Redfield 16 , Caldeira-Legett 17 , and time-convolutionless 
formulations. All these formulations reduce environmen¬ 
tal effects to few terms in £ that contain only the sys¬ 
tem DOF, hence, excluding environmental DOF from the 
explicit consideration. This reduction employs perturba¬ 
tive consideration of system-environment interaction and 
thus is adequate only for weak system-environment cou¬ 
plings. 

Practical necessity to have small system-bath cou¬ 
plings usually requires to consider a large number of sys¬ 
tem DOF in the QME framework, and thus, it increases 
the computational cost of simulation. Therefore, use of 
a time-dependent basis guided by the TDVP for simulat¬ 
ing system dynamics in the QME approach seems as an 
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ideal choice for modelling large open systems. Combining 
the TDVP with QME has been attempted in the past, 
8,18,19 however, these attempts revealed that resulting 
dynamical equations exhibit unphysical behavior which 
was traced to basis set incompleteness. The first prob¬ 
lem is non-conservation of the DM trace, Tr{/S(f)} ^ 0, 
for open systems where only the energy but not the mat¬ 
ter exchange is allowed. 18 ’ 20 The second problem is vio¬ 
lation of energy conservation for an isolated system in a 
mixed state, 6,20,21 Tr{/5(t)£T} ^ 0. It is important to note 
that in the context of the McLachlan TDVP for wave- 
functions, to conserve the energy, there is a requirement 
for a wave-function to be analytic with respect to varia¬ 
tional parameters. This requirement makes the McLach¬ 
lan TDVP 6 equivalent to a principle of least action 22 that 
guarantees the energy conservation. 23 However, violation 
of the energy conservation in the density formalism has 
a different origin and takes place even if the density pa¬ 
rameterization is analytic. It appears as a result of the 
DM entering quadratic.ally in the TDVP, which makes 
Tr{p(t) 2 H} a conserved quantity 20 rather than the en¬ 
ergy E = Tr{p(t)H}. 

These two problems create artificial channels of popu¬ 
lation and energy flows that can lead to inadequate re¬ 
sults. A straightforward solution involving correspond¬ 
ing constraints imposed through the Lagrange multiplier 
method 21 does not preserve a unitary character of quan¬ 
tum dynamics in the isolated system limit. 24 This can 
lead to more subtle artifacts such as non-conservation 
the density matrix purity (Tr{/5 2 }). Another, so-called 
linear mean-field approach has been suggested by Raab 
et al. 19 within the MCTDH framework to resolve these 
issues. However, this method becomes non-variational 
for open systems. 

The aim of this paper is to resolve the non-conservation 
problems using a hybrid approach: To address the en¬ 
ergy non-conservation problem in isolated systems we 
apply the TDVP not to the QME equation but to its re¬ 
cently developed equivalent, non-stochastic open system 
Schrodinger equation (NOSSE) 25 . NOSSE does not have 
issues with energy conservation within the TDVP frame¬ 
work because it is formulated with respect to a square 
root of the density. The root square dynamics is unitary 
for an isolated system and allows for the full reconstruc¬ 
tion of the system density matrix evolution by evaluat¬ 
ing the square at each time. The trace non-conservation 
occurs only for open systems in TDVP NOSSE formula¬ 
tion and is resolved by imposing conservation of the trace 
of the density matrix through a Lagrangian multiplier. 
This constraint does not interfere with the unitarity of 
dynamics in the isolated system limit and thus does not 
introduce any artifacts. 

The rest of the paper is organized as follows. Section 
II is devoted to exposing the non-conservation problems 
at the formal level and details our solutions. Section III 
provides numerical examples illustrating performance of 
our TDVP on a model system. Section IV concludes this 
paper and gives an overview of future work. 


II. THEORY 

A. TDVP problems in the density matrix formalism 


We start by illustrating formally origins of the TDVP 
problems in the density matrix formalism. We introduce 
a finite linear subspace 5? of the total Hilbert space that 
contains time-dependent, not necessarily orthonormal, 
basis functions {\pk(t));k = l,...,Nb}. A projector 
on 5? is built as pN b (t) = I <Pk(t)) [S'~ 1 (t)] fcz |, 

where S k i(t) = (pk(t)\pi(t)) is an element of the over¬ 
lap matrix. The time dependency of the \pk(t)) is given 
through sets of parameters z a k(t ), where the Greek sub¬ 
script labels different parameters corresponding to the 
same basis function. Our only requirement on the pa¬ 
rameterization of the basis is the analyticity condition 22 


d\Vk{t)) = 


( 3 ) 


so that the parameterization preserves the energy in the 
limiting case of an isolated system in a pure state. 

In the time-dependent basis, p{t) is expressed as 


p( t ) = ^2\<Pk(t))B k i(t){(pi{t)\, (4) 

kl 


where Bi k (t) are time-dependent, Hermitian [Bki(t) = 
Bik{t)*\ coefficients. For notational simplicity, when it 
is not essential, the time argument will be skipped in 
the further consideration. A straightforward extension of 
the McLachlan TDVP to the density formalism involves 
minimization of the Frobenius norm of the error 

ii 'p - m\\ = Trw - mpCp - cm 1 ' 2 , w 

which is equivalent to satisfying a stationary condition 20 

Tr{A/5 (3 ,C[p])} =0. (6) 

Replacing p by Eq. (4) in Eq. (6) we obtain EOM for B = 
{Bki} and z a k (detailed derivation is given in App. V) 


B = S~ 1 LS ~ 1 - (. S-WB + Bt^S- 1 ) , 

S/3 1 Ckl z pl = P k > 

where [r\ki = (p k \tpi), [L]ki = {p k \B[p]\ pi }, and 

dpi 


_ / dpk 

kl \dz ak 


1 -Pi 


N b 


1-Pi 


dzip 


[.BSB ] 


Ik, 


N b 


( 7 ) 

( 8 ) 


( 9 ) 


C[p\ pi ) Bik■ (10) 


Equation (8) is solved for i = {z ak } as a linear equation 
where z and Y are vectors and C is a matrix 


z = C~ l Y. (11) 

A given parameterization can lead to redundancies be¬ 
tween some of the parameters z ak and B k i (e.g., in the 
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MCTDH and vMCG methods). In such cases extra rela¬ 
tions between redundant parameters must be combined 
with Eqs. (7) and (11) to eliminate any redundancies. 

Equations (7) and (11) are equivalent to EOM ob¬ 
tained in the particular case of moving Gaussians de¬ 
rived in Ref. 8. This is not fortuitous since the “type II” 
density parameterization of Ref. 8 coincides with that of 
Eq. (4). By choosing the MCTDH ansatz for the time- 
dependent basis, it can be verified that Eqs. (7)-(ll) pro¬ 
vide the EOM obtained in Refs. 19 and 20 (up to the 
removal of redundancies in parameterization). 

To illustrate variation of the total system population 
we will use the derived EOM to consider time derivative 
of the system population 

Tr{,3(t)} = Tr{S5 + (r + r t )B}. (12) 

Using Eq. (7) the variation of the DM trace on time is 
given by 


Tr{/5(t)} = Tr {S^L - S^rBS - + (r + t^)B) 

= Tr {S-'L} = Tr {P Nb C[p}}, (13) 

where in the second equality we have used cyclic invari¬ 
ance and linearity of the trace operator. For incom¬ 
plete bases Ppf b C[p] ^ £[/5], and even if Tr{£[p]} = 0, 
Tr{P/v b £[/5]} can be non-zero for open systems as has 
been shown by Raab et al. for the Liouvillian in the 
Lindblad form. 20 In contrast, for isolated systems the 
population is always conserved because 

Tr {P Nb C[p}} = —i Tr{S~ 1 (HBS - SBH)} 

= -iTr{[H,B}}=0, (14) 

where [H] k i = {ip k \H \pi). 

For an isolated system, the energy variation with time 
is given by 

Tr{j5H} = Tr {BH + BH} 

= Tr {HB - BHS^t - S" 1 HB} ,(15) 


where in the second equality we used the B definition 
from Eq. (7), [H] k i = {<p k \H\<pt) + {ip k \H\<p t ), and cyclic 
invariance of trace. Equation (15) can be further simpli¬ 
fied using hermiticity of B, H , and S 


Tr{m = J2 2 Re 

akl 



dipk 

dz ak 


= -2Im[Tr{i t Y'}], 


r ^ 

A 

\ 

i - p N b 

H 

) Bi k 


(16) 


where 


[Y] 


a 

k 



i-P Nb 



Bik- 


(17) 


Using Eq. (11) and the vector notation for i and Y we 
arrive at 


(18) 


Equation (10) simplifies for isolated system to 



i -Pn„ 



[BSB]i k . (19) 


According to Eqs. (16) and (18), energy conservation 
takes place at least in four special cases: 


1. Time-independent bases, where z = 0. 

2. Complete bases, where Pjv 6 = 1, and therefore, 

Y = 0. 

3. Systems in pure states, where BSB = B (the 
idempotency condition), then Y = Y and us¬ 
ing hermiticity of C, it is easy to see that 
Im[Y'tC , - 1 Y] = Im[ytC'- 1 Y] = 0." 

4. Linear parameterization, \<p k ) = )T) n ,\<j> a ) z ak , 

where \4> a ) are time-independent basis functions. 
In this case, Eq. ( 8 ) becomes CqzBSB = 
Y 0 BSB , where [C 0 ] Q/ 3 = (</>a|[l - Pv b ] \4>p) 

and [Y 0 ] a i = - PN b ]H\<pi), which 

leads to i = C^Yo, and similarly Y = 
Y q B. Equation (16) gives — 2 Im[Tr{idY}] = 
— 2 Im[Tr{l^Cb' 1 lQ-B}] = 0 , where we used her¬ 
miticity of Cq and B , and semi-positivity of B. 

We would like to conclude this section by emphasizing 
that a general non-linear time-dependent parameteriza¬ 
tion (e.g., used in MCTDH and vMCG methods) of the 
density matrix will not conserve energy of the isolated 
system in a mixed state. 


B. Constrained TDVP for NOSSE 

To resolve the non-conservation problems we will 
use a problem specific approach: The population 
non-conservation for open system will be corrected 
by introducing a constraint using Lagrange multiplier 
method. For the problem of energy non-conservation 
we reformulate dynamical equations using non-stochastic 
Schrodinger equation for a square root of the density. 
This difference in approaches to the two problems stems 
from an objective to preserve unitarity of the dynamics in 
the isolated system limit. Imposing a constraint ensuring 
energy conservation breaks the unitarity of the isolated 
system dynamics and may lead to artifacts. 24 

To proceed to formulation of EOM for a density square 
root in open systems, we introduce some notation first. 
Matrix B is semi-positive, and therefore, there exist a 
matrix M such that B = M. Defining 

I m k ) = ^ I Vi) M ik , 

1 


Tr {’pH} = -2Im[l rt C'- 1 U] 


( 20 ) 
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and using Eq. (4), we can rewrite the density matrix as 
P = M M o-kM£ k (<p b | 

abk 

= ^2\m k ) {m k \. ( 21 ) 

k 

All states (|mfc)} can be arranged in a vector form m = 
(\rrii) \ 1 ri 2 ) . ..), such that m represents a density square 
root because p = mm). In Ref. 25 we found EOM for 
m for a system with a given Liouvillian operator C 

ra = 7C[m], (22) 

where /C[m] is a vector of states that satisfies a Sylvester 
equation 

K[m\m t + m/C[m]^ = C[mm)]. (23) 

Note that for an isolated system 7C[m] reduces to —iHm 
and Eq. (22) becomes an uncoupled set of TDSE for each 
component of m. This reduction guarantees that ap¬ 
plication of TDVP to Eq. (22) will be free from energy 
conservation issues arising in the density formalism. In 
order to optimize the parameters of m, the Frobenius 
norm ||m — 7C[m]|| is minimized. 

To address non-conservation of the total population 
for open systems we introduce a corresponding constraint 
using the Lagrange multiplier method. The Lagrangian 
function is 

A = ||m — /C[m]|| 2 + aJ^ Trjmm*}, (24) 

where A is a Lagrange multiplier. Minimizing A with 
respect to parameters M = { M k {\ , z, and A leads to a 
constrained TDVP 


where we employed the Sylvester equation [Eq. (23)] in 
the matrix form, L = KM^ S + SMK t. 

Now, it is easy to show that Eqs. (28-31) con¬ 
serve all required quantities. Taking the trace of p = 
■§i Sfc/ Wk) B ki (tpi | and using Eq. (31) gives 

^Tt{p} = Tt{BS+(t + t^B} 

= Ty{S- 1 L-Ti{S~ 1 L}BS} =0, (32) 

where the last equality holds since Tr{/5} = TrjfJS 1 } = 1. 
Therefore, the trace is conserved for any Liouvillian L[p\. 
If the system is isolated, C[p\ = —i[H,p] in Eq. (31), it 
is possible to show that a trace of any power of p is also 
conserved 

^ Tr^} = Tt{BS(BS )"" 1 + B{t + r t )( J B5)”“ 1 } 

= -iTrlS-^HBS - SBHjtBS)"- 1 

- Tr {S~\HBS - SBH)}(BS) n } = 0. (33) 

This illustrates unitary character of dynamics. The en¬ 
ergy is also conserved for the isolated system limit 

4 Tr{ pH} = Ty{BH + BH} 
ot 

= Tt{BH - S^tBH - HBt ] S'- 1 } 

= -2 Im [z)Y] = -2 Im [Y^C^Y] = 0, (34) 

where we invoke the hermiticity of C to see that 
Y^C~ 1 Y does not have imaginary component. 

III. NUMERICAL EXAMPLES 


Re [Tr { 5rh t (m — 7C[m] + Am)}] = 0, (25) 

J^Tr{mmA} = 0. (26) 

Replacing m by its explicit form [Eq. (20)] in Eq. (25) 
and Eq. (26) we obtain EOM for M and z (details of the 
derivation are given in App. VI) 


M = V'Jf-V’rM 

— i Tr{KM^ + MK^}M (27) 

i = C~ 1 Y, (28) 


where K ki = {ip k \Ki[m]) and 
/ dip k 


r< a P — , _ 

kl \dz a 


y k = 


E ! dp k 

\dz ak 


1 -Pn, 
i -Pn„ 


< 291 


\[lCi[m]) M* kl . (30) 


Using the relation B = MM t + MM\ Eq. (27) can 
be recast into an EOM for B 


B = S~ 1 LS~ 1 - (S _1 tB + Bt^S- 1 ) - Tr{S _1 T}B, 

(31) 


A. Model 


To illustrate performance of our developments for the 
case where system quantum effects are affected by in¬ 
teraction with environment, we will consider the sim¬ 
plest linear vibronic coupling (LVC) model 26 of crossing 
surfaces. 2 '’ 28 The system contains two electronic states, 
donor | D) and acceptor |A), which are coupled through 
two nuclear DOF 


^ = Ey(?« + ^) \\ D ) W + W < A \ 


Q! = l 


(35) 


dir [ |£>> (D\ - | A) {A\ ]+cx 2 [ | D) (A\ + \A) (D\ 


where oj a is the frequency of the coordinate x a , and d 
and c are diabatic coupling constants. Numerical values 
of the parameters are taken from Ref. 29 for the two- 
dimensional model of the bis (methylene) adamantyl rad¬ 
ical cation: w 1 = 7.743 • 10~ 3 a.u., u >2 = 6.680 • 10~ 3 
a.u., d = 5.289 • 10 -3 a.u., c = 9.901 • 10 -4 a.u. The dis¬ 
sipative part is introduced through bilinear coupling of 
the system coordinates x a with coordinates of harmonic 
bath oscillators. System-bath couplings are taken small 
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enough so that the second order perturbation theory can 
be applied for the system-bath interaction in conjunction 
with the rotating wave and Markovian approximations. 
The resulting Liouvillian 30,31 has a Lindblad form 15 

2 

C[p\ = - i[H , p] + hiY (2 Z a pZl - ZlZ a p - pZlZ Q ) 

a=l 

2 

+h 2 Y (?ZlpZ a ~ Z a Zip - pZ a Zi) , (36) 

a—1 

where hi and h 2 are parameters that depend on the bath 
temperature, and Z a operators are defined as 

Zi = (“■' |D> <D|+ ( Sl+ ^71) |3> {Al 

Z 2 = ai(\D)(D\ + \A)(A\), (37) 

with annihilation operators a a = ( x a + ip a )/y^2. For 
C[p\ in the Lindblad form, K[m\ is derived in Ref. 25, 
and for our model is given by 

2 

K.[m] = —iHm + hi ^ ^ ZapZ^m J) -1 — Z^Zam'j 

Oi—1 

2 

+h 2 Y [ZaPZairn))- 1 - Z a Zlrnj . (38) 

a=l 

B. Numerical details 

We apply our constrained TDVP to the varia¬ 
tional multiconfiguration Gaussian (vMCG) ansatz. 10 In 
vMCG, the time-dependent basis consist of products 

I <Pk) = Wk ) | gk) (39) 

of discrete electronic parts \ak) = \D) ,\ A) and 

spatial two-dimensional moving Gaussian parts | gk) 
parametrized as follows 

(x\g k ) =exp(-x T A k x + ^x) , (40) 

where a; is a vector of the mass- and frequency-weighted 
nuclear coordinates, Ak are frozen widths (Ak = 0) and 
£k are parameters encoding the position (Re [£*,]) and mo¬ 
mentum (Im[£fc]) of a Gaussian. This choice of parame¬ 
terization satisfy the analyticity condition d | ifk) /d£ k a = 
0 . 

For reasons of numerical stability it is more convenient 
to work with coherent states (CSs) of harmonic oscillators 
rather than frozen Gaussians. 11,32 However, EOM de¬ 
rived from the McLachlan TDVP expressions for CSs can 
have energy conservation issues because CSs do not sat¬ 
isfy the analyticity condition. On the other hand, TDVP 
EOM for frozen Gaussians can be projected onto the ba¬ 
sis of CSs, {| Zk)}, of harmonic oscillators with frequencies 
ui a by multiplying Eq. (40) with Vfefc = e~^ k Re ^ fc l/ 2 / 7r 
and setting Ak = I 2 / 2 . The transformation of Gaussians 


| Zk) = Vkk | gk) leads to a transformation of the matrix 
B = VBV\ where [V] M = V kk hi- In the CS repre¬ 
sentation variational parameters become z ak = Cfca/v^, 
they correspond to eigenvalues of the annihilation oper¬ 
ator a a for CSs \z k ) 

a a | z k ) = z ak | z k ) ■ (41) 

Upon the V transformation, the structure of the 
Eqs. (28-31) remains the same, thus, we will use these 
equations to propagate parameters of CSs. 

The population loss for an open system is given by 
Eq. (13). In the particular case of the Liouvillian defined 
by Eq. (36) in the basis of CSs, the expression for the 
population loss becomes 

2 

TrM = 2 Y Tr{[hiZl(P Nb - 1 )Z a (42) 

a=1 +h 2 Z a (P Nb - l)Zl]p}. 

Equation (42) can be further simplified using the defi¬ 
nitions of Z a [Eq. (37)], p [Eq. (4)], Pjv b , and applying 
Eq. (41) 

2 

Tr{p} = 2h 2 Y Tr,: V {da(Av f) - 1 )cP a 

“ =1 x[(D\p\D) + (A\p\A)]},m 

where Trjv{...} is the trace over system nuclear DOF. 
Thus, the hi containing term disappears and the popu¬ 
lation loss takes place only if h 2 ^ 0 . 

The initial system DM p is taken as a Boltzmann dis¬ 
tribution pb in the donor electronic state 

-P D HP D /k B T p 

Pb =- . . . ; -, (44) 

Tr { e -PDH p D /k B T p D y 

where Pd = \D) (D\ is the projector on the donor elec¬ 
tronic state, ks is the Boltzmann constant and T is the 
temperature. An approximate initial density matrix p(0) 
is obtained by minimizing 

A = ||m — m.B || 2 + A[Tr{mm^} — 1], (45) 

where m b is a square root of pb■ This minimization 
is equivalent to a two-step procedure: first, maximizing 
Tr-fpB-PiVi,} with respect to the parameters of the basis, 
and second, constructing the initial density in the opti¬ 
mized basis as p( 0 ) = P Nb p B pN b /Tr{p B P Nb }- 

Solving the EOM (28-31) requires inversion of the 
S and C matrices, these matrices can become singu¬ 
lar and lead to numerical difficulties . 33 This problem 
is addressed by applying a regularization procedure 34 
on singular eigenvalues (s) before the matrix inversion: 
s —> s + ee _s / £ , where e = 10 -6 is the chosen threshold. 
This procedure provides a faster and more stable propa¬ 
gation with a numerical precision of the order of e. The 
limiting step of our implementation is the inversion of the 
C matrix. The computational cost of this inversion scales 
cubically with the dimensionality of C. For the system 
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with d DOF described by Nt, Gaussians, the dimension¬ 
ality of C is 2Nbd (where the factor of two accounts for 
position and momentum parameters of Gaussians), and 
thus the computational cost scales as 8 N^d 3 . 

Simulations with time-dependent basis are compared 
with exact simulations that are done by projecting 
the QME onto the direct-product basis of the two- 
dimensional harmonic oscillator. 360 static basis func¬ 
tions are employed to obtain converged results that we 
will denote as p e (t). 


C. Isolated system simulations 

For an isolated system where hi = h 2 = 0, C[p\ = 
—i[H,p], and K[m] = — iHm , we choose the tempera¬ 
ture of the initial Boltzmann distribution to be T = 1000 
K (or T = 0.47wi = DA1uj 2 ). Figure 1 shows population 
dynamics of the donor state, Tr{pPjj}, simulated with 
the NOSSE-based formalism. A good agreement with 
the exact propagation is obtained already for 25 CSs, 
and results of our method converges to those of the exact 
propagation with 35 CSs. Root-mean-square deviations 
(RMSDs) from the exact dynamics for various properties 
are defined as 

RMSD(A) = ^J q * dtTi{A[p(t)-p e (t)]}\ (46) 

where tf = 100 fs, and A = 1, Pd, and H for the total 
population (Pt), the donor state population (Pd), and 
the system energy (E), respectively. Table I shows that 
the system energy in QME-based simulations is not well 
conserved. As a result, the corresponding donor popu¬ 
lation dynamics have larger deviations than those in the 
NOSSE formalism. 

The coupling between two electronic states is relatively 
weak and the population dynamics can be analyzed using 
perturbation theory. Fast scale oscillations are produced 
from transitions between (m,n)D and (m, n± 1 )a levels 
due to the linear coupling [cx 2 in Eq. (35)] , here, the first 
and the second numbers correspond to the numbers of vi¬ 
brational quanta along the tuning ( x \) and coupling ( x 2 ) 
modes, and the subscripts denote the electronic states. 
Vibrational levels (m, n)u and (m, nil)^ are separated 

TABLE I. RMSDs of the donor population and energy with 
the NOSSE and QME approaches for different number Nb of 
CSs. 


N b 

15 

25 

35 


RMSD(E), 

in units of ui 


NOSSE 

< 10" 5 

< itr 5 

< 10 -5 

QME 

2.8 • 1CT 3 

2.3 • 10 -3 

2.4 • 1CT 3 


RMSD (P D ) 


NOSSE 

5.9 • 1(T 3 

1.1 • 10 -3 

1.7- 10" 4 

QME 

6.1 • 10 -3 

1.1 • 10 -3 

4.9 • 10" 4 



FIG. 1. Isolated system: Population dynamics of the donor 
electronic state Tr{p(t)pD} for the exact and vMCG dynam¬ 
ics based on constrained NOSSE [Eqs. (29) and (30)] formal¬ 
ism with different number Nt, of CSs. 


by the energy difference ±w 2 , which translates into the 
population oscillation period t p = 2n/u} 2 = 23 fs. There 
are also slower time-scale oscillations that appears as a re¬ 
sult of transitions between the (m, n )d and (m=Fl, nil) a 
levels. The tuning mode does not contribute to the cou¬ 
pling between electronic states, but it shifts minima of 
these states so that levels’ couplings are modulated by 
the Franck-Condon factors. The energy difference be¬ 
tween the (m,n)D and (m=Fl,n±l)^ levels is ±(uq — 1 x 2 ) 
that results in 143 fs time difference between two subse¬ 
quent maxima of slow population oscillations. 


D. Open system simulations 

We simulate the system (T = 1000 K) located in the 
donor well and interacting with a harmonic bath accord¬ 
ing to Eq. (36) with hi = 0 and h 2 = 3.675 • 10 -4 a.u. 
For /i 2 / 0 the trace of p varies if no constraints are 
applied [see Eq. (43)]. In the constrained NOSSE-based 
vMCG simulations, the system density trace is conserved 
up to 2 • 10 -5 . The population dynamics of the donor 
state Tr{p(t)PD} is given in Fig. 2, it illustrates that 
the constrained NOSSE vMCG formalism gives the ex¬ 
act dynamics for 35 CSs. A comparison of deviations 
from the exact dynamics in trace-constrained and un¬ 
constrained NOSSE simulations is given in Table II. As 
expected, the unconstrained NOSSE donor populations 
have larger deviations than those of the constrained for¬ 
malism. Although generally there is no notable compu¬ 
tational cost difference between constrained and uncon¬ 
strained schemes, in some cases the constrained scheme 
has better efficiency due to more regularly behaving so¬ 
lutions of the corresponding EOM. Thus, all open system 
calculations similar to the ones reported here should be 
done with the constrained formalism. 

The fast-scale oscillations (23 fs period) observed in the 
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FIG. 2. Open system: Population dynamics of the donor elec¬ 
tronic state Tr{p(t)PD} for the exact and vMCG dynamics 
based on constrained NOSSE [Eqs. (29) and (30)] formalism 
with different number N b of CSs. 

isolated system case are also present in Fig. 2, and their 
nature is exactly the same. However, due to influence of 
the environment, the long term dynamics has changed in 
the open system case. The origin of the irreversible decay 
of the donor population is a two-step-resonance process: 
( 0 , 0 )d —► ( 0 , 1)d —> ( 0 , 0 )^ 4 , where the first step is driven 
by the environment excitation of the system within the 
donor state, while the second step is due to the linear 
coupling between the electronic states. Since there is no 
energy difference between the initial ( 0 , 0)d and the final 
( 0 , 0 )a states, this process leads to the larger population 
transfer observed in Fig. 2. 


problem, the developed approach can significantly re¬ 
duce the number of propagated parameters for simu¬ 
lating non-adiabatic dynamics while conserving all re¬ 
quired physical quantities. Moreover, even larger effi¬ 
ciency improvements compare to standard techniques are 
expected in molecular problems with a large number of 
nuclear DOF. Although the current implementation of 
the constrained TDVP-NOSSE approach has been illus¬ 
trated using vMCG ansatz, working equations can be 
easily used to derive other system density parametriza- 
tions (e.g., MCTDH). For a better numerical stability we 
have used a basis of CSs that are localized in space. Thus 
our implementation is perfectly suited for a combination 
with the on-the-fly evaluation of electronic potential en¬ 
ergy surfaces . 32,35,36 From the application prospective, 
our developments open a new venue of the on-the-fly 
photochemistry with incoherent light. The sun light is 
incoherent, hence, realistic modelling of chemical photo¬ 
stability, solar energy harvesting and utilization requires 
introducing light as the system environment in a mixed 
state, which can be easily handled by the current ap¬ 
proach. 
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IV. CONCLUDING REMARKS 

In this work we exposed the energy and total popu¬ 
lation non-conservation problems that occur in the den¬ 
sity matrix TDVP formalism. The proposed constrained 
TDVP-NOSSE approach resolves the non-conservation 
issues by combining the Lagrange multiplier method to 
preserve the total population and the NOSSE formal¬ 
ism to reproduce unitary dynamics in the isolated sys¬ 
tem limit. As illustrated on a model surface-crossing 

TABLE II. RMSDs of the total and donor populations with 
the constrained and unconstrained NOSSE approaches for dif¬ 
ferent number Nb of CSs. 


V. APPENDIX: DERIVATION OF THE EQUATIONS OF 
MOTION FOR THE TDVP APPLIED TO QME 


Main steps in derivation of Eqs. (7) - (11) are detailed 
below. Considering the parameterization of the DM in 
Eq. (4), variation of Sp in Eq. ( 6 ) can be split into con¬ 
tributions from Bki, z a k , and z* fe 


S P=J2 


kl 


dp 

dB k i 


SB k i 


r dp 


y , — 

y ^dz ak 

ka. 


Sz n 


dp 


■Szl, 


(47) 


N b 

15 

25 

35 

Constrained 

RMSD(Pr) 

< 1CT 5 < 1CT 5 

< 10" 5 

Unconstrained 

8.8- 10“ 2 

2.2 • 10” 2 

6.5 • 1CT 3 

Constrained 

RMSD (P D ) 

1.3 • 10“ 2 1.9 • 10~ 3 

6.2 • 10 -4 

Unconstrained 

6.8- 10“ 2 

1.8 • 1(T 2 

5.0 ■ 10" 3 


Owing to independent character of SBu, Sz ak , and Sz ka 
variations in Eq. (47), Eq. ( 6 ) represents a system of 
equations 


Trj 

[asl (p r[^])| 0 , 

(48) 

Trj 

[&(3-£0)}= 0 , 

(49) 

Tr ■ 

[*(V^])} = 0 . 

(50) 

















( 54 ) 


Replacing p by its explicit form [Eq. (4)] into Eq. (48) 
leads to 

E { E ( I Vc) B cd {ip d | + | ip c ) B cd (tp d \ 

ab cd / i \ 1 

+ I<A=) B cd (<p d \J -£[p]| I <Pa) 

= ^ ) {^SkcBcd.S d l “h TkcBcdBdl T SkcBcdTld^j P kl — 0? 

cd 

(51) 


where Sk c = (^fcl^c) is the overlap matrix between time- 
dependent basis functions, Tk c = {Vk\<Pch and B ki = 
(ipk |£[/5]| ipi). Using the matrix notation, Eq. (51) is com¬ 
pletely equivalent to Eq. (7). 

Equations (49) and (50) are complex conjugates of each 
other, and therefore, constitute one unique stationary 
condition. Replacing p by its explicit form [Eq. (4)] into 
Eq. (49) leads to 


E S ° b 


ab 


dtp b 


dz<y.k 


{£( | <Pc) B cd (ip d I + I ip c ) B cd (if d | 

kc) B cd {ip d \) - £[/3]j I ip a ) = 0. 


cd 


note that variations 


v dM u 


and 


8m) = 


\ - dm) 


4^ 3M* SMkl 


kl 


\ - dm 1 * 

^dz*J ak 


(55) 


kcx 


span the same tangent space of m(t) because varia¬ 
tions of parameters ( 8M kl , 8z* k ) and their time deriva¬ 
tives ( 8M kl , <5i* fe ) are completely arbitrary. Therefore, 
Eq. (25) can be substituted as 

Re [Tr {5ml (m — 7C[m] + Am)}] = 0. (56) 

The parameterization of m(t ) is chosen so that the basis 
functions | tpk) are analytic with respect to their param¬ 
eters, d\ifik) tdz* ak = 0. This analyticity allows us to 
extend the zero condition in Eq. (56) from a real part to 
the total trace expression 22 

Tr {(5ml (^ _ JC[m] + Am) } = 0. (57) 


Substituting B cd by Eq. (7) we obtain 
/ dip k 


J2\~gT L , {j2( K \<P c )BcdS d a + \<Pc)B cd (tp d \ip a ) 

a ' cd . 

+ \<p c ) [S-'LS- 1 - S~ 1 tB - B^S-\ d Sda) 

- C[p\ I ifia) | Bak 

= 53 ( | { E ^ ^ I Vc) BcdSda) 

a ' 1 cd . * 'i 

+ (-P/V b — l)i3[p] | ip a) jB a k 


= o, 


(52) 


The variation <5mA is split into contribution from Mki and 
z a k [Eq. (55)], and using independence and arbitrariness 
of variations for different parameters, Eq. (57) is rewrit¬ 
ten as a set of equations 



Replacing ml by its explicit form [Eq. (20)] in Eq. (58) 
leads to 


where we used the definition of the projector p N b = 
Du \vk) [s- _1 ]fcj (ipi |. The time derivative in Eq. (52) can 
be expanded using the chain rule in terms of the basis set 
parameters d/dt = Dip Zipd/dzip , thus 


/ difk 

^\dT 


ak 


i -Pi 


N b 


£) [BSB ^ 


/ dVk 
\ dz, 


ak 


i p 


N b 


P[p] Wi) [ B]ik 


= = o, 

i/3 


(53) 


dM* ( / 

E QM* { E ( l^ c ) Mco, + |<p c ) M ca 

ab kl c \ 'i 

+A | <p c ) M ca J - |/C a [m]) | 
= E ( 'SkcM c i + TkcM c i + A SkcM c ^j - Ku = 0, (60) 

C 

or using a matrix notation for Eq. (60) gives 

M = S~ 1 K-S~ 1 tM-\M. (61) 

Inserting Eq. (61) into the constraint Eq. (26) gives the 
expression of the Lagrange multiplier 


that is Eq. (8) with definitions given in Eqs. (9) and (10). 


VI. APPENDIX: DERIVATION OF THE EQUATIONS 
OF MOTION FOR THE CONSTRAINED TDVP APPLIED 
TO NOSSE 

Here we provide main steps in derivation of Eqs. (27) - 
30). Starting from the constrained TDVP Eq. (25) we 


N Tr {Mfft + JfMl} 

' “ 2Tr{MtSM} 

= i Tr {MK^ + KM^}, (62) 

where the second equality arises because the 
Tr{ M^SM} = Tr{mmt} = Tr{/3} = 1. Substi¬ 
tuting A by Eq. (62) in Eq. (61) leads to the final 
equation of motion for M [Eq. (27)]. 
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To derive the EOM for z a k we replace mA by its explicit 
form [Eq. (20)] in Eq. (59) 




+A I ip c ) Mca'j - \JC a [m}} | 


= 0. 


(63) 


Substituting M ca by Eq. (61) gives 
/ dip k 


J2 M *ka 


\dz, 


E l 9(fk 
\dz a k 


{ ]T ( Wc) [S~ l K - S~ 1 tM - A M} c 
<Pc) M c a + A | ip c ) M c 2j - I IC a [m}} | 

{^(i-P Nb )\<j> c )M ca 
C +(P Nb - 1) I K, a [m]) )M* ka 


= 0. 


(64) 


Expanding the time derivative using the chain rule in 
terms of the basis set parameters gives 


/ 9(pk 

W' dz<xk 


l-Pi 


N b 




j dtpk 


\ dz, 


ak 


i -Pi 


N b 


\1Ci[m])M kl 


= E<^%- y *“ = 0, 

W 


(65) 


where we used the definitions from Eqs. (29) and (30). 
Adopting the matrix notation the last equation can be 
rewritten as Eq. (28). 
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